{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5000000/5000000 [00:28<00:00, 175726.93it/s]\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "class langevin:\n",
    "    def __init__(self) -> None:\n",
    "        self.dN = 1.0\n",
    "        self.dt = 0.01\n",
    "        self.dt2 = self.dt * 0.5\n",
    "        self.dt4 = self.dt2 * 0.5\n",
    "        self.kT = 1\n",
    "        self.x = 0\n",
    "        self.v = 1\n",
    "        self.mass = 1.0\n",
    "        self.k_spring = 1\n",
    "        self.f = -self.k_spring * self.x\n",
    "        self.Ek2 = 0.0\n",
    "\n",
    "        self.temperature_coupling = 100\n",
    "        self.c1 = np.exp(-0.5/self.temperature_coupling)\n",
    "        self.c2 = np.sqrt((1 - self.c1**2) * self.kT)\n",
    "        self.c2m = self.c2 * np.sqrt(1/self.mass)\n",
    "\n",
    "    def langevin(self):\n",
    "        self.v = self.c1*self.v + self.c2m*np.random.normal()\n",
    "    \n",
    "    def run(self, nsteps, sample_intervals):\n",
    "        self.length = int(nsteps/sample_intervals)\n",
    "        self.x_list = np.zeros(self.length, dtype=float)\n",
    "        self.v_list = np.zeros(self.length, dtype=float)\n",
    "\n",
    "        for step in tqdm(range(nsteps)):\n",
    "            # half step\n",
    "            self.langevin()\n",
    "            self.v += self.dt2 * (self.f / self.mass)\n",
    "            self.x += self.dt * self.v\n",
    "            # another half step\n",
    "            self.f = -self.k_spring * self.x\n",
    "            self.v += self.dt2 * (self.f / self.mass)\n",
    "            self.langevin()\n",
    "\n",
    "            if step % sample_intervals == 0:\n",
    "                nn = int(step/sample_intervals)\n",
    "                self.x_list[nn] = self.x\n",
    "                self.v_list[nn] = self.v\n",
    "\n",
    "lan = langevin()\n",
    "lan.run(5000000, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAokAAADUCAYAAADwWLkhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7K0lEQVR4nO3dd3hUZfr/8fedSW8kQALSA4QuKIQmiID0ahd7XcS17Or6tazrrnV1f8uqa9lVxIJlRVYsIFgIIgqI0lEIUUSEUEJCkZKQNs/vjwxsSJkJJDPPlPvlNRfJnDnnfAzcmWdOuR8xxqCUUkoppVRFYbYDKKWUUkop/6ODRKWUUkopVYUOEpVSSimlVBU6SFRKKaWUUlXoIFEppZRSSlWhg0SllFJKKVWF9UGiiDhEZI2IfGQ7i1KBQkRGiUi2iGwWkXurWT5YRH4VkbWux59t5FQqUGmNKQXhtgMAvwOygETbQZQKBCLiAJ4HhgM5wAoRmWOM2VjppV8ZY8b5PKBSAU5rTKlyVo8kikgLYCww3WYOpQJMH2CzMWaLMaYYmAlMtJxJqWCiNaYU9k83Pw3cDTgt51AqkDQHtlf4Psf1XGX9RWSdiHwsIl19E02poKA1phQWTzeLyDhgjzFmlYgMdvO6ycBkgLi4uF4dO3byTUCl3Fi9elW+MSbF0u6lmucqz6+5GmhtjDksImOAD4D0ajdWqcY6ddIaU/atWhUcNab1pfxRbevL5jWJA4AJruKKBhJF5E1jzJUVX2SMmQZMA+jVK8Ms/Wal75MqVUlMhPxicfc5QMsK37cAdlZ8gTHmYIWv54vIv0SksTEmv/LGKtZYRkaGWblSa0zZJxIcNab1pfxRbevL2ulmY8x9xpgWxpg2wCTg88oDRKVUtVYA6SKSJiKRlNfPnIovEJGmIiKur/tQXut7fZ5UqcCkNaYU/nF3s1LqJBhjSkXkVuBTwAG8YozZICJTXMtfAC4CbhaRUqAQmGSMqXy6TClVDa0xpcr5xSDRGPMF8IXlGEoFDGPMfGB+pedeqPD1c8Bzvs6lVLDQGlPK/t3NSimllFLKD+kgUSmllFJKVaGDRKWUUkopVYUOEpVSSimlVBU6SFRKKaWUUlXoIFEppZRSSlWhg0SllFJKKVWFDhKVUkoppVQVOkhUSimllFJV6CBRKaWUUkpV4RfT8oW6us72aTj1DQhSt5172r53N6+UUkopL9EjiUoppZRSqgodJCqllFJKqSp0kKiUUkopparQQaJSSimllKpCB4lKKaWUUqoKHSQqpZRSSqkqdJColFJKKaWq0D6J9cRTr0N3vQw9rluHbXvmft0wD40OPbdBPPVGidpjUSmllLJHjyQqpZRSSqkqdJColFJKKaWq0EGiUkoppZSqQgeJSimllFKqCh0kKqWUUkqpKnSQqFQAEpFRIpItIptF5F43r+stImUicpEv8ykV6LTGlNJBolIBR0QcwPPAaKALcJmIdKnhdX8DPvVtQqUCm9aYUuW0T2Itee5l6P4FTjeLy9wtBEqdzlPeNrjPJh76GDrC3C8Pd7hf7m71ME8fUYyHHo2h20exD7DZGLMFQERmAhOBjZVedxswG+jt23ihq7S0lOLiYmJjY21HUXWjNeanCgsLiYiIIDxchy++oEcSlQo8zYHtFb7PcT13nIg0B84HXvC0MRGZLCIrRWRlXl5evQYNJR1GXUtkVDRxcXE0aNmBjBsfY/yzS44/VECptxrT+qofRUVFdOvWjdjYWBITE5k0aRKbN2+2HSvo6SBRqcBT3THUyoeLnwbuMcaUedqYMWaaMSbDGJORkpJSH/lCUnLa6aQNvpgOo6/HWVrCqpf/xIbZz3g8y6D8Ur3VmNZX/YiKimLixIk8+uijXHfddcybN49evXqxcOFC29GCmh6vVSrw5AAtK3zfAthZ6TUZwEwpPyffGBgjIqXGmA98kjBEzJs3j88//5ypU6eS2qUvqV36ApA+8mqyPvw38U1aIyF8XUQA0xrzA8YYbrrpJq644gqmrndA09HlCxKh711D+f7dp3n0i1ye3lh+pH7ubQMtpg1OOkhUKvCsANJFJA3YAUwCLq/4AmNM2rGvReQ14CN986pf27Zt46qrrqJNmzYUFRWdsEzCHHQ5/9bj3ztLS3wdT9WN1pgfePHFF3nppZdIS0uD+LNPWBbTsAm9Jz9+/HutMe/Q081KBRhjTClwK+V3VGYBs4wxG0RkiohMsZsuNBhjuOGGGygpKeGdd94hOjq6xtfu/m4JXzx2Bfv37/dhQlUXWmP2bdmyhTvuuIORI0dyzz331Pg6Ywzr3nqc796Z6sN0ocPakUQRaQm8DjQFnMA0Y8w/beVRKpAYY+YD8ys9V+0F9MaYa32RKZR88MEHZGZm8uyzz5Kenu72tTHJTSjYl8tDDz3E008/7ZuAqs60xuy66667cDgcvPzyy4S5aYUhIkTEN2DLwrdZsWIFvXvrjeb1yeaRxFLgD8aYzkA/4Jbq+lAppZQ/cTqd/N///R/dunVjyhTPB5UatEin1Vnjee6558jOzvZBQqUC29KlS3n//fe5//77ad68ucfXp4+4hqiEhvz+97/XG8XqmbUjicaYXcAu19eHRCSL8hYDlftQ+SiPp+XuX1DmYXlJac3Lj5a4vwG1qNR9n8SSMvfL3fVh9NQHMdLh/nNEdITD7fKoiJrXD/fUo9Hj9f7aR1H5XlhYGHPnzqWwsLDWvdo6jrmB/DUL+Otf/8qMGTO8nFCpwNavXz/ee+89Ro8eXavXR8TE0WH0dSyb9Q8WLlzIsGHDvJwwdPjFjSsi0gY4E/immmWTgckALVu18m0wpZSCmvscLq1d/8OohGRuuukmnnnmGR5//HGaNWtWj+mUCnxVayyFV15aWev1W/Qdw94lbzN16lQdJNYj6zeuiEg85R3rf2+MOVh5+Qk9phprjymllD2713/Jypf/RPGRX0963bvvvpuvvvpKB4hKubHuP0+wZdE7J72eIyKSWbNm8frrr3shVeiyeiRRRCIoHyC+ZYx5z2YWpZTyZMuiWRw9kEdETPxJr9u0aVOaNm3qhVRKBYeCvTvZ/s180kdcfUrrDxgwoJ4TKWtHEqW8A+nLQJYx5klbOZRSqjYO7d7Kvp/W0XrgRCTM/bW4NSksLOTGG2/U6xKVqsa2rz8ChFZnTTjlbXz77bcMGzaMffv21V+wEGbzdPMA4CpgqIisdT3GWMyjlFI12rZsLuIIp0Wf2l1MX53o6GhWrlzJU089pXdhKlWBs6yU7cvnkdqlHzHJqae8nejoaBYuXMhrr71Wf+FCmLVBojFmiTFGjDHdjTFnuB7zPa+plFK+VVZSRM6KT2l6+tlEJSSf8nZEhClTprBu3TpWrFhRjwmVCmx7Niyj6OA+WtfhKCJA9+7dOeuss5g2bZp+EKsHfnF3sz8wVeZuP5GbLjKA+xY3AAXFpTUuO3S05mUAhwvdL99fVOx2ebGbFjlRDvenzRpERbhdnhjjfnlCTM3/xGIj3e9bPPTACfPwd+auRY62x1Enw1lSTJuB55HSqU+dtzVp0iR+97vf8dZbb9GnT923p1QwiG6QQuuB55Pimv+8Lq699lomT57M6tWr6dWrVz2kC13W725WSil/FxGbQMexN9KwXfc6byspKYnx48czc+ZMSkvdfwBUKlQkte7M6ZfcSZij7seuLrroIiIjI3nzzTfrIVlo0yOJSinlxsGDB8n9fikpnfoQFu7+yLk7FfvA7UntT2wnGPdk5gl3Ss+9bWCdsioViFauXMnBnVtIbNa2XraXnJzMnXfe6XHKTOWZDhKVUsqNDz74gBXT7mXAHf8mOa1bvWwztXNfUjvX/bSaUsHgvvvuY+XaLIY88DZSh2uBTmjI3Wws3x+B9ys16dYPYidHTzcrpZQbM2fOJKZhU5LadK3X7RpnGft+Wk9ZiftripUKZrm5uXz++ec07zWsTgPE6pQUHmb/L1Zm+g0aOkhUSqka5Ofns2DBApr1PLfe38Dysley7J+3kP/DqnrdrlKB5N1338XpdNKs57n1vu3vZj3JihfvwTjL6n3boUIHiUopVYO5c+dSWlrKaWcMqfdtN2p/JuFRsexet7jet61UoJg9ezZdunQh4bS0et9209MHUnz4APu2fFfv2w4VOkhUSqkaLF68mBYtWtCgZYd637YjIpLUrv3J/X6pHulQIamgoIAVK1YwceJEr2w/tUs/wsIj2b3+S69sPxSEzI0rHntqelhe5qFRYlGp+1/y7noh5h4scrvuT78edrs8O6/Q7fK8QzVf89Q4PtLtuh1Sot0u75ic4Ha5IarGZWEezt6FeTi9F+ahj6JSdfXKK6+wY8cOfvvBL17ZftMe57Bz9UL2bfmORu3P8Mo+lPJXsbGx7Nq1i6KiIq79T1a9bz88OpbGHTPYve5Lupx/W71fMhIK9EiiUkrVICwsjJYtW3pt+6ld+hIWHknu90u9tg+l/Fl8fDyNGjXy2vZP6zGIwv25HNr5k9f2Ecx0kKiUUtV44IEHuPvuu726j/CoWAbc8S86jZvs1f0o5W+Ki4sZOnQo8+bN8+p+mvY4h3P++AaJzdt7dT/BSgeJSlkkIrNFZKyIaC36EafTycsvv8zPP//s9X01aNmxTk26lXtaY/7pyy+/ZNGiRZSVefd63IiYeBKatvHqPoKZFo1Sdv0buBz4UUSeEJFOtgMpWLVqFbt27WLChAle35ezrJSsOS+wY+UCr+8rRGmN+aE5c+YQHR3NsGHDvL6vgzu3sOb1Ryg6uM/r+wo2OkhUyiJjTKYx5gqgJ7AVWCAiy0TkOhHRw0uWfPTRR4SFhTFmzBiv7yvMEU7u90vZ/s18r+8rFGmN+R9jDHPnzmXYsGHExsZ6f39lJexY+Rl5m771+r6CjQ4SlbJMRBoB1wI3AmuAf1L+hqaHlixZsGABvXv39uoF9RWldunH3s1rOXzYfScDdWq0xvzLTz/9xNatWxk9erRP9pfYPJ2oxIbs2fi1T/YXTHSQqJRFIvIe8BUQC4w3xkwwxrxjjLkNiHez3igRyRaRzSJybzXLJ4rIehFZKyIrRUQnLK0lYwzdu3dn0qRJPttnapf+mLJSFi5c6LN9hgqtMf9TUFDAmDFjGD58uE/2J2FhpHbuR17Wt5SW1tyOTlUVMn0SPfHQBpFSp9Pt8qMl7pcfdtMncetB90cPFm8+4HZ51s/ur7PYu7egxmVxce7PtuzqmOp2eVgH932nYiIcNS6LdrMMICrc/V9KRHD0SZxujDnhPKOIRBljiowxGdWtICIO4HlgOJADrBCROcaYipOULgTmGGOMiHQHZgF6LVYtiAgvvPCCT/fZsO3phEfFMn/+fK81Fg5hWmN+pnv37l6/q7mylC792P7NfJYvX87AgTqery09kqiUXY9W85yncyJ9gM3GmC3GmGJgJnDCyMIYc9iY4y3k4/DYLl4dk5ubi/HYfb9+hYVHcNoZg4mMdN/cXp0SrTE/UlZWxp49e3y+35ROvYlv0poDBw74fN+BTAeJSlkgIk1FpBcQIyJnikhP12Mw5afF3GkObK/wfY7rucr7OF9ENgHzgOvrJ3nwGzJkiE9PNR/T44r7ePbZZ32+32ClNeaf1q5dS5MmTZg7d65P9xsRE8/g+99k3LhxPt1voNPTzUrZMZLyC+lbAE9WeP4Q8EcP61Z3nr3KUQxjzPvA+yIyCHgEqLbXhIhMBiYDtGrVylPuoLZjxw6ysrK4/np77/dHjx4lOtr9dJiqVvyixrS+TpSZmQlAnz59rOy/rKwMp9NJRITe2F4bOkhUygJjzAxghohcaIyZfZKr5wAV54prAex0s68vRaSdiDQ2xuRXs3waMA0gIyMjJE+ZjX92CQA5334CwJy8Rix2PedLY8eORUT46KOPfL7vYOMvNab19b/6Alj+6rskNGvLjbN+BH70aY7De7aTkjKBF154gUsuucSn+w5UerpZKQtE5ErXl21E5M7KDw+rrwDSRSRNRCKBScCcSttvL67Z7EWkJxAJ7K3n/42gk//DKiLjGpDYrJ2V/bdu3ZovvviCkpISK/sPJlpj/qesuIh9P62ncYdq7xfyuthGp+F0Oo8fzVSe6SBRKTviXH/GAwnVPGpkjCkFbgU+BbKAWcaYDSIyRUSmuF52IfC9iKyl/C7NS42v78YIMMYY8rJX0rhjLyTMzq/G4cOHc+TIEZYvX25l/0FGa8zP7P/5O5ylxTTuaGeQGOYIZ8iQISxYsMDnN6cFKj3drJQFxpgXXX8+dIrrzwfmV3ruhQpf/w34W10yhhxj6H7JXUTEJVqLMGTIEMLCwsjMzOTss8+2liMYaI35n4TT2nL6pf9Ho/Y9rGUYNmwYH3zwAVu2bKFdOztnDAJJyAwSjYfuBJ4+VZR5aKRYUuq+T+KBozWfPvpud6HbdX/cfsDt8p821XipDACFebtrXhjj/g0xKsr9P5EmDdxfYJ+WGFfjsgYx7i8c9tSbssy475MYLv7bR1FEnnG33Bhzu6+yqHISFkaT0wdYzZCUlETv3r1ZsGABDz10SmMb5aI15n+iEhvSeoD350N351gD7wULFuggsRZCZpColJ9ZZTuAOtHO1QuJS21FgxbpVnPcddddHD161GqGIKE15kdKCg6xa92XND19AJHxSdZypKen89BDD9G3b19rGQKJDhKVssB156XyE86yUta//f9oljGc7pfeZTXLRRddZHX/wUJrzL/k/7ia9W8/QXyT52locZAoIvz5z3+2tv9Ao4NEpSwQkaeNMb8XkblU33/N7jmZEPPrtk2UFhXQuEMv21EA+PHHH9m9e7del1gHWmP+JT97JY6oGJJad7EdhZKSEpYtW0abNm1o3bq17Th+TQeJStnxhuvPqVZTKADysleCCI079LQdBYDbb7+drVu3kpWVZTtKINMa8yN52Stp1P5Mwhz2hx379+9n8ODBPPbYY/zxj576qoc2bYGjlAXGmFWuPxdTPo/sfmAf8LXrOeVD+dkradCiA5FxDWxHAcovrt+0aRM5OTm2owQsrTH/UbB3FwV5OdZa31SWmppKjx49tF9iLeggUSmLRGQs8BPwDPAcsFlERttNFVqKi4s5uGOz35xqhvI2HVB+B6aqG60x+w78Un5EPMVPBolQ/kFs6dKlFBQU2I7i1+wf9w0SZR5a6JS4aedy6Gip23X37XPfIqdw51a3y9nnpkVOUhO3qx461NLt8oLiMrfLj5bVvNzTz8xD1yGP3LU9kmqnZrXiH8AQY8xmABFpB8wDPraaKoRERkYy/NEPKSspsh3luNNPP53U1FQWLlzIddddZztOoNMas6xZz6E0bNedqMRGtqMcN2zYMKZOncqSJUsYMWKE7Th+y+qRRBEZJSLZIrJZRO61mUUpS/Yce/Ny2QLssRUmVDkio4i02ES7MhFh6NChLFq0SGeGqDutMT8Q3aAx4ke9awcOHEhERASff/657Sh+zdqRRBFxUD6V0XDKJ1NfISJzjDEbbWVSyldE5ALXlxtEZD4wi/I7MC+mfN5Y5SNXXXUV201LWvYdYzsK459dcvzrgo4X0K3HlUx4bunx5+beNtBGrICkNeYfNmzYwIpp99F54hTim9i/k7hijfW/cxrfNW1zwnNaYyeq1ZFEEYl2TYr+nojMFpE7RMT9VBue9QE2G2O2GGOKgZnAxDpuU6lAMd71iAZygXOAwUAekGwvVmjJzc3lzTffpOjXvbajVBHb6DS/uZEmQGmN+YFPP/2U3O+X4IiMsR2lisTm7f3ibmt/VtufzuvAIeBZ1/eXUd5e4OI67Ls5sL3C9zlAlRboIjIZmAzQslWrOuxOKf9hjNELzfzAsVNN/nLXZWXbl8+jcP8eOozWfy4nS2vMP2RmZhLfpDUxyam2o1RRUniYHz95jdQu/fz2d4BttR0kdjTGVJyRe5GIrKvjvqu7OKG6hqfTgGkAvXpl6MU5Kqi4jsjfAHSl/IgHAMaY662FCiGZmZkkJSXRoGUH21GqtX/rRnauyqT9iKv0iMcp0hqzp7i4mMWLF5OaMcp2lGo5IqPZtmwuZcVHdZBYg9reuLJGRPod+0ZE+gJL3by+NnKAirfOtgDc3IarVFB6A2gKjAQWU14Hh6wmChHGGBYsWMDQoUORMIftONVq3DGD0qICft22yXaUQKY1Zsny5cspKCigcQf/HICFOcJp2P4M8n7Qab5rUttBYl9gmYhsFZGtlDcmPUdEvhOR9ae47xVAuoikiUgkMAmYc4rbUipQtTfGPAAccc01OxY43XKmkHDkyBHOOOMMxo8fbztKjRqnnwlAvr6J1YXWmCVFRUX06dOHRu3PsB2lRo07ZlCQl0PBvt22o/il2p6/qPdjxcaYUhG5FfgUcACvGGM21Pd+6ounvnoOD7f2x4TXfKSiSUKk23WTk93fI7SjUTO3y8vcHSVJaOx23cTEKLfLE6LdH4GJdtS83NPPzNNyT38nftQL0Z0S158HRKQbsBtoYy9O6IiPj2fOnPLPpbMr3N3oTyLjk0hskU5+9irSR15jO06g0hqzZPjw4QwfPvyEu4f9zbEm+vnZq2jVf6zlNP6nVoNEY8wv3ti5MWY+MN8b21YqQEwTkWTgAcqPpMe7vlZedvDgQRIT/ac3Yk1Su/TjwC9Z2i/x1GmNWVBcXIzT6SQ6uq6NULwr4bQ0Ek5rS1mx+0krQpVeCa2URcaY6a4vFwNtbWYJJWVlZaSlpXHzzTfz6KOP2o7jVqdxk21HCGhaY3Z8+umnXHLJJXz99de2o7glIpxz3wzbMfyWzt2slEUi0khEnhWR1SKySkSeFhH/mbsqSK1Zs4Z9+/bRtWtX21FqzTjdT4Gpqqc1ZkdmZiYiQufOnW1HqTWtsap0kKiUXTMpnyLsQuAiIB94x2qiEJCZmQnA0KFDLSepne9n/5OlT/3WdoxApTVmQWZmJoMGDSIqyv117f6gpPAInz9yGT8vnm07it/RQaJSdjU0xjxijPnZ9XgUSLIdKthlZmbSvXt3mjRpYjtKrUTGNeDAtizy8/NtRwlEWmM+tnPnTjZu3MiwYcNsR6mViJg4BO0iUB0dJCpl1yIRmSQiYa7HJcA826GCWWFhIUuWLOHcc8+1HaXWGnfoBcawaNEi21ECkdaYjy1cuBAgYAaJUF5jezevoaSkxPOLQ4gOEpWyQEQOichB4CbgP0Cx6zETuMNmtmBnjOH555/nyiuvtB2l1pJadyY8Kvb4m6/yTGvMnr59+/LEE0/QvXt321FqrXHHDMqKClmxYoXtKH4lZO5u9thTz0NLPU/LI8Pdj7cTIiNqXNYq2X2fxC5tGrpd7qkzxp7clBqXNUhy356gQ6skt8vTG7tfPyGy5n9inn5mYR4+wnj6OzFVZ3n837qWeygaYxKsBghhsbGx3HDDDbZjnJRjM0Mcu5ZSeaY1Zk+HDh245557bMc4KY3Se4IImZmZnHXWWbbj+A09kqiUZSIyQUSmuh7jarnOKBHJFpHNInJvNcuvEJH1rscyEelR3XZC0ezZs9m+fbvtGCetVf9xTJkyhbIyvQPzZGmN+c6OHTuYO3cuBQUFtqOclMi4RNJHXkP//v1tR/ErOkhUyiIReQL4HbDR9fid6zl36ziA54HRQBfgMhHpUullPwPnGGO6A48A0+o7eyDau3cvF198Ma+++qrtKCetafezueuuu3C4mcVIVaU15luzZ89mwoQJ5Obm2o5y0jqOuYHhw4fbjuFXQuZ0s1J+agxwhjHGCSAiM4A1QJUjFxX0ATYbY7a41pkJTKT8DRAAY8yyCq9fDrSo59wBaeHChRhjGDFihO0op2Tv3r1s3ryZvn372o4SSLTGfOizzz6jffv2pKWl2Y5y0owxbNy4kdjYWNq0aWM7jl/QQaJS9iUB+1xfN6jF65sDFc+X5gDuRg03AB+fUrIgUHHe2HVvv0F4TDyPLD9K2Ar/nU+2JrfffjsLFy5k165diKeLclVFSWiNeV1xcTFffPEF11wTmPOMO0uKOfPMM7ntttuYOnWq7Th+QU83K2XXX4E1IvKa6wjHKtdz7lQ3Oqj2Lh0RGUL5G1iNV5GLyGQRWSkiK/Py8moZO/AYY8jftILGHXoR5gjMz8fDhg0jNzeXDRs22I4SSKzWWKjUF8DXX3/NkSNHAvaUrSMyigEDBmgXgQoC8zelUkFARMIAJ9AP6E35G9M9xpjdHlbNAVpW+L4FsLOa7XcHpgOjjTF7a9qYMWYaruupMjIyPNwrH7gK9+2mcH8u7YcHTuubyo71dszMzKRbt26W0/g/f6ixYK+vikfqf/jkDSTMwQs/RPPys4F3pB7Ka+xPf/oTeXl5pKTU3BkkVOiRRKUscV0jdasxZpcxZo4x5sNavHkBrADSRSRNRCKBScCcii8QkVbAe8BVxpgf6j18AIptdBrDH/2QZr0Cp8FvZa1atSI9PV1b4dSS1phvpY+4mnP++AYRMfG2o5yyYw3AtXF9OT2S6OLp8p4Ih/vxdHSk+zsO46JrXt6lYaLbdcM8hGvuoddh7qHiGpfFR7nP3bVprNvlHZPdtyJLiKm5P6Snn1m4h0aJnq7Ist0LsZYWiMhdlM8le+TYk8aYfTWtYIwpFZFbgU8BB/CKMWaDiExxLX8B+DPQCPiX69q1UmNMhvf+NwJDVKL7nqOB4Nxzz+XNN9+kpKSEiIia60sdpzXmIxIWRnxqS88v9GO9evUiMTGRzMxMLrnkEttxrNNBolJ2XU/5tU6/rfR8W3crGWPmA/MrPfdCha9vBG6sp4wBz1lWyurXHqTN2eeXT3EXwO666y7uuOMOwsP113ctaY35wJ6sb9i1ZhGdJ/6WyDj3Bz78WXh4OJ988gmdOnWyHcUv6OlmpezqQnk/tnXAWuBZoKvNQMHowNaN7F63mJLCw7aj1Fm7du3o0KGD3t1ce1pjPrB73WJ2rf2C8Gj3Z58CQf/+/UlOTrYdwy/oIFEpu2YAnYFnKH/z6ux6TtWjvE3fgoTROL2n7Sj1Yv78+Tz66KO2YwQKrTEvM8aQt2kFjTr0DNjOARUVFxfzxBNP8PHHId/VSAeJSlnW0RhzozFmkesxGehoO1SwycteSVLrzkTEBsd0vosXL+bhhx/myJEjnl+stMa87EheDoX7dpPSsbftKPUiIiKCZ555hhkz9LOEDhKVsmuNiPQ79o2I9AWWWswTdIoLDnHglyxSOgXHGxjA8OHDKSkp0Tswa0drzMvys1cCBE2NiQjDhw9nwYIFIT9Xug4SlbKrL7BMRLaKyFbga+AcEflORNbbjRYcig/tI6lVJ1I69bEdpd6cffbZxMXFMX/+fM8vVlpjXiYiNGzXg9jGzW1HqTdjxoxh3759fPvtt7ajWBX4Fw/UkqdrvD1dAu6hGwsRDvdbSIg+9VYV3R3uZ5FqEV9zixuAMlNz/9bwMPe5k6Mi3S531+KmfHnN/8Siwt3/UB0esnn6SwuQ6/pH2Q4Q7OKbtGbgH160HaNeRUVFMWzYMObNm4cxRm9icU9rzMtaDzyP1gPPsx2jXo0YMQKHw8H8+fPp37+/7TjWhMwgUSl/ZIz5xXaGYOZ0OiktKiA8KvDvuKxs7NixZGdnk5+frzNDuKE15l2lRQU4IqIRT0dSAkxycjKDBw9m374a22mGhOD6W1VKqQpWrVrFZ/eOI891zVQwueGGG8jKytIBorJq09xpLHrkMowz+K7d++yzz3j++edtx7BKjyQqpYLWvHnzcDpLadC8ve0o9WJ8NfPhGqfzhKM4c28b6MtIKoQZY9iz4Wvim7ZGwtzPoBUotMZOpEcSlVJBa968eSS37kpkfJLtKF6xY1UmC+6fQHHBIdtRVAjKzs6mYO9OUrsE7zV7q179C6tf+4vtGNboIFEpFZR2797NypUrSe0avG9gMcmpFB/5lfxNK2xHUSFo3rx5ADTpepblJN4TEZvAnqxvcJaW2I5ihQ4SlVJB6dhsCcE8SExu05WI2ET2bPzadhQVgubNm0fCaW2JadjEdhSvadKlP2VFhez7aZ3tKFboIFEpFZQGDRrEP/7xDxKD5HrE6kiYg5TOfdizcTnG6bQdR4WYe+65h45jb7Qdw6sadehJWHgkezYutx3FCr1xxUU8NN3z0AaRSA89/9xt31M7wOgI9xcEN4h136vQWXObRMI9/I95Wu4pm7teiJ56NHrqqODp70yFtnbt2nHnnXeyqJoL0YNJk65nsXNVJgd+2UhyWjfbcVQIGTlyJM/9EGc7hleFR8XQKP1Mdn+/lM7n3RJyPUl1kKiUCjpr1qzhl19+YezYsbajeF1ql360HTqJyPhk21FUCJk5cyZdunSxHcMn2px9PgV7d2GcZYgjtIZNofV/q5QKCc8//zz//e9/ycvLsx3F6yJiE+hy3i22Y6gQcvToUX7zm99w+eWXQ7erbMfxuibdBtiOYI2VaxJF5O8isklE1ovI+yKSZCOHUir4lJaW8uGHHzJu3DgiI91PKxksnGWl5P+wioK9u2xHUSEgMzOTw4cPc/7559uO4jPFRw6ya+0XtmP4nK0bVxYA3Ywx3YEfgPss5VBKBZklS5aQn5/PBRdcYDuKz5QUHGL583eyffk821FUCHjvvfdITExk6NChtqP4TM6KT1n1ygMc3rPddhSfsjJINMZ8ZowpdX27HGhhI4dSKvi89957REdHM2rUKNtRfCYqIZlG7Xqwe/2XtqOoIFdaWsqcOXMYP358yBypBzit+yAAdq8LrRrzhxY41wMf17RQRCaLyEoRWZmXH/zXFyml6mbNmjWMHDmSuLjgvuuysqY9BnFo189kZ2fbjqKC2KZNmygsLAypU80AMQ2b0KBVZ3avX2w7ik95bZAoIpki8n01j4kVXnM/UAq8VdN2jDHTjDEZxpiMlMY6kb1Syr0vv/ySGTNm2I7hc01dRzref/99y0lUMOvWrRv5+fmMHz/edhSfO63HIA78ksX27aFzytlrdzcbY4a5Wy4i1wDjgHONMW46+fmG59ZHHnr6eVg7MtxNn0QPK0dFuH9BmbtGiIC7n26Yh16FDg8/GIen9d0sr2sfxBBrV6VqwRiDiNCgQQPbUXwuJjmVpNZdWLhwIffee6/tOCoIHauvmJgY21GsaNpjEJvmvsiiRYu4+uqrbcfxCVt3N48C7gEmGGMKbGRQSgWX0tJSunbtyksvvWQ7ijW9rn/4+HSEStW3jz/+mB49evDTTz/ZjmJFfGorhv7lnZAZIIK9PonPAVHAAlf38uXGmCmWsiilAtT4CrOp7Mn6hqysLF78Jo85R4N7lpWaxCQ3ITxc298q73jrrbfIycmhZcuWtqNYE9uome0IPmXlt4kxJngnU1XKB1xH4/8JOIDpxpgnKi3vBLwK9ATuN8ZM9X1K39q5KpPwmHhSu/SzHcWqGTNmMH36dBYvXkyYp2s6VI20xsod+yBWWlTIgnffo3nvkVz44reWU9ljnE4uv/xyunTpwp/+9CfbcbxOf4MoFWBExAE8D4wGugCXiUjl+bH2AbcDQfnGVVlZ8VF2rVvMaT3OwREROm05qhMeHs6SJUtYunSp7SgBS2usqtzvl1JWfJTmvdzebhD0JCyMffv28fLLL+N0Om3H8TodJCoVePoAm40xW4wxxcBMYGLFFxhj9hhjVgAlNgL6Wu73SykrKqR5xgjbUaw777zziIuL480337QdJZBpjVWyY+VnRCen0rBtd9tRrLvyyivZunUry5Ytsx3F63SQqFTgaQ5U7MGQ43rulJzQizRA5zqOb5pG2yGX0qh9D9tRrIuLi+P8889n1qxZHD161HacQFVvNRYM9QXQPGMEHUZdj+glDJx33nnExsbyxhtv2I7idfq3rVTgqa75zym3kTqhF2lKYPYiTWzWli7n34qEOWxH8QvXXnstBw4c4N1337UdJVDVW40FQ30BNO81jFb9x9qO4Rfi4+O5+OKLefvttzl8+LDtOF6lt8HVkqeefGEeevoZqfn3S4S4H6uHexjKG0+/u+rQhVI8/I97aJPotr2k9kE8ZTlAxdsLWwA7LWWxbvf6r4hq0Ijk1pUvGQtdQ4YM4ZZbbqFDhw62owQqrTEX43Sy9av3aNbzXKISkm3H8Ru33HILzZs3p6QkuK820EGiUoFnBZAuImnADmAScLndSHY4y0r57r9P0qB5e/pM+bvtOH4jLCyM5557znaMQKY15pKfvZINs/9JZHxSyN+0UlHv3r3p3bu37Rhep6eblQowxphS4FbgUyALmGWM2SAiU0RkCoCINBWRHOBO4E8ikiMiifZSe0fexm8o+jWfVmeF3hRhtZGdnc2HH35oO0bA0Rr7n21fzyUiNvH4tI/qf8rKyvjoo4/YtGmT7Sheo0cSlQpAxpj5wPxKz71Q4evdlJ8iC2o/L/4v0Q1SSO16lu0ofqNig/HVr/2FvKxvOffh9wiPKp9Kbe5tA21FCyhaY7Bt2zZ2r/+KtMEXh3xrqeocOnSISy65hMsvv5zp06fbjuMVOkhUSgWk9evXk//DKjqNn0KYQ3+VVafNoAvZufpztn8zn7RBF9qOowLMs88+C0DaoIssJ/EvFT+IpfQcwasz3mBb2gSiEhsefz5YPozp6WalVEDatGkTUQ0a02rABNtR/FZy2ukktenKz4tmYZxltuOoALNz505OO2MwMQ2b2I7it9oOvgRnaTFbl7xvO4pX6CBRKRWQLrnkEs598L9ExibYjuK3RIS2Qy6lYO9Odn8XmvNZq1P31ltvccZVwT/1XF3EN2lFk24D+OWr9ykrDr6+pDpIVEoFnE2bNmGM0dPMtXBaj0HEN21DQX5IdnBRp6C4uJjNmzcDaI3VQtuhk3BExXAkf4ftKPVO//briaeefu56AnpsY+hh23Vog+iRp16GHtfXXoeqnuXn55ORkcFtt90GzbS5rycS5mDQPa/qm72qtRkzZnDzzTezevVq21ECQsN2PRjywNtBWWN6JFEpFVCmTp1KQUEBV199te0oAePYm9eBbeVHYJWqSXFxMY8++igZGRmcfvrptuMEBBEhzBFOWUkRB3dsth2nXukgUSkVMPLy8njuuee47LLL6Ny5s+04ASX3u6Usmfob5s+f7/nFKmS99tprbNu2jQcffNDjjFvqRGvfeIxv/n0XZcVFtqPUGx0kKqUCxsMPP0xhYSEPPPCA7SgBJ6VLX2IbNeOBBx7A6XTajqP80OHDh3nwwQfp168fI0eOtB0n4LQ++3yKDu4NqjuddZColAoIR48e5dNPP2XKlCl06tTJdpyAE+YIp+O437BmzRpeeeUV23GUH1q2bBn79u3jySef1KOIp6Bx+pmkdO7Lj5+8Rm5uru049SL4rrJUSgWNik1rAdpP+Tc/lZZWeV7VTrOe5xK/5XPuu+8+LrzwQpKTk21HUn5kxIgRbN++nZSUFNtRAlbXC25n8RPXcN999wXFhzEdJCql/N6v238gvkkrHJHROCKibMcJWCLCM888w4gRI9iwYQMDBwbHrBDq1Bz7sGWMYf/P39GwbXfXkmx7oQJcfJNWpA2+mOzsbIqKioiKCuzfV3q62Q+I1PHh4b8wOfVHXbMpVVdFh/bzzb//wNo3H7MdJSicccYZbNu2TQeI6rgdKxew7OlbyP1uqe0oQaHjmBv56quvAn6ACDpIVEr5uQ2zn6ak8DDpo661HSVoREdH43Q6mTZtGocOHbIdR1l09OBeNsx+muQ23Ujt2s92nKDgiIgkLCyMPXv28Prrr9uOUyc6SFRK+a2cbz9h5+rP6TDqOhKbtbMdJ6isW7eOKVOmcMcdd9iOoiwxzjLWvflXykqK6HHFfUiYw3akoPL3v/+da6+9lkWLFtmOcsp0kKiU8kvr1q1j/TtTadj+DNoNu9x2nKBz5plncu+99/Lyyy8zffp023GUBT98/Cp5m76l6wW3E9+kle04Qecvf/kL6enpXHrppWzfvt12nFOiN64opfxSUlISKR17033S/wXldFe2VLwz3DQZSUqnTCZPuZnXNpaSnNYNgLm36fWKoSChWTvaDLqQVmdNsB0lKMXHx/PBBx/Qt29fzjvvPJYsWUJMTIztWCdFjyQqpfzKkSNHKC0tpXXr1vSe/DhRiQ1tRwpaEubgzGsfIia5Cate+TNlJcEzU4Sq2a+//gpAszOH0O2i32tPRC8Y/+wSxj+7hLsz99Jp0v2sXrOGrmOuOf58oNCP50opv1FQUMCECRNo2LAhs2bNsh0nJETGJtD35qkU7N2l7YVCwMqVKxk5ciQvvfQSkGo7TkhocvoAel7zII079rId5aTpkUSllF/Iz8/n3HPPZdGiRUycOFGPbvhQXEoLUjr1BspvFlq9erXlRMobPvvsM4YMGUJCQgJ9+vSxHSekNOs5lMi4BjhLS/jh41c5cuSI7Ui1ooPEIKC9ClWg27BhAwMGDGDNmjXMnj2bK6+80nakkFRWUsQPH7/KoEGD+OCDD2zHUfXEGMOLL77I2LFjadu2LcuWLaNFixa2Y4WkvT+t44dPXmPw4MHk5OTYjuORnm5WSllx7LocZ1kpix65nLKSo2RMeZJXclJ4JYCu2Qkmjogozvr98+z/4FHOP/98fvvb3zJ16tSAu9henXiDUv4Pq1n+3O9I6dSH0654iJtmbwG22AsXwlI6ZtD7xr+y8T+P0L17d6ZPn84FF1xgO1aN9EiiUsqKQ7t+xllaQpgjnJ7X/IVB97xKw3bdPa+ovCq6QWOWLl3KnXfeyb/+9S969uxJcXGx7VjqJBmnk4M7NgPQKP1MMn7zOH2m/J2ImHjLyVST0wewevVq2rVrx4UXXsiDDz5oO1KNdJColPKptWvXcumll7L4iWvY+uVsAJLTuhKd2MhyMnVMVFQU//jHP/j888+56aabiIyMBGDRokU4nU7L6ZQ7paWlvPXWW3z5t2tZ8uRNFO7LRURoevpAJEzf8v1Fhw4dWLZsGY899hgTJpS3INq+fTs//PCD5WQnsnq6WUTuAv4OpBhj8m1mUSqQiMgo4J+AA5hujHmi0nJxLR8DFADXGmOs3Y1w7Jqo119/na+//pqEhATanXs5LfqOthVJufG/U5UR4Mhg4bNL2P/z9yx96mZiGzenecZw5jx1Dx07dgzaG4wCrcZ27tzJU089xX/+8x927txJwmlp9Lj8PqKTUmxFUm4cr7GEQXy9tACWLuG7d6byy7I5NO6QQfOMYXzx3N0kJSVZzWltkCgiLYHhwDZbGZQKRCLiAJ6nvH5ygBUiMscYs7HCy0YD6a5HX+Dfrj+9rqysjM2bN3P5EzMpOrSftHMuAmDZP1+gpOAQnSfeTKv+44mITfBFHFVPGrTqxJnX/IVtX8/lx09n0Lnza7Ru3ZqPP/6Yzp07c+DAAaKjo4mOjrYdtc78vcaOHDnC+vXrWbt2LZ06dWLIkCEUFxfzzDPPMGrUKG688UZe3NJAjxwGmA6jrycyIZkdKz5l3VuP03jm/2PcuHHHbyLbs2cPKSkpPv1gZvNI4lPA3cCHFjMoFYj6AJuNMVsARGQmMBGo+AY2EXjdGGOA5SKSJCKnGWN2uduwMYaCggJKSkqOP1JSUggPDyc/P5+cnBx+/fVXDhw4cPzPKVOmEBkZyVNPPcW0adP45ZdfKCwsBMARGUOrs8bjiIii9+S/ER4dG7RHnoJdmCOc5r2G0bzXMI7+ms/udV+yd/Ma/jB/J47MvWTNeYEtn88kpmFTYpKbMLZ/N1q1asXDDz+MiLBu3Tr2799PYmIiUVFRREZGEhMTc/wu22PXPTocDsLsD268VmPFxcUn1JeI0LhxYwB+/PFHDhw4cEJ9NW3alHHjxgEwfPhwsrOzycnJoXy30KLvaM64IgKAIY98iDMmnmlbQaz/CNXJikpsSMcxN9Bh9PUc2LqRIXE5xy/1AOjZsycHDx4kLS2NVq1a0bJlS0aMGMF5550HwCeffEJiYiKxsbFERUURFRVFw4YNSUpKwhhDUVERDocDh6P2c3RbGSSKyARghzFmnb5hKHXSmgMVJwLNoeoRjOpe0xxw+wa2evVq4uLiTnhuy5YtpKWlMX36dO67774q63x0OI2ohGS2r9zDgagmNO3bncRmbUls0YGEpm0ICy9/A4uIiauyrgpM0Q0a02bQBbQZ9L+7Mpt0PYswRzhH8nIo3J/LrLmfUHq0gLWpIwFYPeNhdq5acMJ2UlNTyc3NBeDiiy9mzpw5AFx00UU++j+pkVdqLCsri6ioExuWDxgwgCVLyk89Tpw4kaysrBOWp3TqTd+fk8rXPyDIaV1J7zacxObtSWyRTkxyk+Ov1ZtSgoOIkJzWlbV0BcpPTRtjaHj25UTs3ELuvt1sXZNN4edfMv+7Xby8vTHOslLm31H18p0//OEPTJ06lUOHDtGgQYOTz3Ls00h9E5FMoGk1i+4H/giMMMb8KiJbgYyarkkUkcnAZNe33YDvvRC3PjUG/Pn6Sn/PB4GRsaMxxsr5UhG5GBhpjLnR9f1VQB9jzG0VXjMPeNwYs8T1/ULgbmPMqmq2F0g1Fgj/NjRj/QiKGguw+oLA+LehGeuuVvXltSOJxphh1T0vIqcDacCxo4gtgNUi0scYs7ua7UwDprnWXWmMyfBW5vrg7xn9PR8ETkaLu88BWlb4vgWw8xReAwRWjfl7PtCM9SVYaiyQ6gs0Y33x94y1rS+fX7VgjPnOGJNqjGljjGlDeaH1rG6AqJSq1gogXUTSRCQSmATMqfSaOcDVUq4f8Kuna6WUUsdpjSmFzriiVMAxxpSKyK3Ap5S353jFGLNBRKa4lr8AzKe8NcdmyttzXGcrr1KBRmtMqXLWB4muo4m1Nc1bOeqRv2f093ygGT0yxsyn/E2q4nMvVPjaALecwqb9/Wfv7/lAM9aXYKwx/bnXD81Yd7XK57UbV5RSSimlVODSTkpKKaWUUqqKgB0kishdImJEpLHtLBWJyN9FZJOIrBeR90UkyXamY0RklIhki8hmEbnXdp7KRKSliCwSkSwR2SAiv7OdqToi4hCRNSLyke0s3uKv9QVaY6cqUOoLtMZs89ca8+f6guCssYAcJIp/T+m3AOhmjOkO/ABU7T5sgfxvmqnRQBfgMhHpYjdVFaXAH4wxnYF+wC1+mBHgd0CWx1cFKD+vL9AaO1WBUl+gNWab39VYANQXBGGNBeQgkf9N6ed3F1QaYz4zxpS6vl1Oee8sf3B8miljTDFwbJopv2GM2WWMWe36+hDl/4Cb2011IhFpAYwFptvO4kV+W1+gNXaqAqG+QGvMH/hpjfl1fUFw1ljADRKlwpR+trPUwvXAx7ZDuNQ0hZRfEpE2wJnAN5ajVPY05b/cnZZzeEWA1RdojZ0SP64v0BrzN/5SYwFTXxA8NWa9BU51pBZT+vk20Ync5TPGfOh6zf2UH3p+y5fZ3Khukmy//BQrIvHAbOD3xpiDtvMcIyLjgD3GmFUiMthynFPm7/UFWmPe5K/1BVpjvhSANRYQ9QXBVWN+OUisryn9fJ3vGBG5BhgHnGv8p8dQradps0lEIigvrreMMe/ZzlPJAGCCiIwBooFEEXnTGHOl5Vwnxd/rC7TGvMXP6wu0xnwmAGvM7+sLgq/GArpPoohsBTKMMX4zibaIjAKeBM4xxuTZznOMiIRTfgHyucAOyqedutwYs8FqsAqk/LfmDGCfMeb3luO45foEdpcxZpzlKF7jj/UFWmOnKpDqC7TGbPLHGvP3+oLgrLGAuyYxADwHJAALRGStiLzgaQVfcF2EfGyaqSxglj8Vl8sA4CpgqOtnt9b1aUepirTGTo3Wl6otv6uxAKgvCMIaC+gjiUoppZRSyjv0SKJSSimllKpCB4lKKaWUUqoKHSQqpZRSSqkqdJColFJKKaWq0EGiUkoppZSqQgeJSimllFKqCh0kKqWUUkqpKnSQGCJEpLeIrBeRaBGJE5ENItLNdi6lgoHWl1LeJyJtRGSTiMxw1du7IhJrO1cw02baIUREHqV8rsYYIMcY87jlSEoFDa0vpbxLRNoAPwMDjTFLReQVYKMxZqrdZMFLB4khREQiKZ/v8ihwljGmzHIkpYKG1pdS3uUaJH5pjGnl+n4ocLsx5jybuYKZnm4OLQ2BeMrn5Iy2nEWpYKP1pZT3VT6ypUe6vEgHiaFlGvAA8BbwN8tZlAo2Wl9KeV8rEenv+voyYInNMMEu3HYA5RsicjVQaoz5j4g4gGUiMtQY87ntbEoFOq0vpXwmC7hGRF4EfgT+bTlPUNNrEpVSSinl91zXJH5kjNHOAT6ip5uVUkoppVQVeiRRKaWUUkpVoUcSlVJKKaVUFTpIVEoppZRSVeggUSmllFJKVaGDRKWUUkopVYUOEpVSSimlVBU6SFRKKaWUUlX8f67jCEsBD3QtAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 792x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3)\n",
    "axs[0].hist2d(lan.x_list, lan.v_list, bins=21, range=[[-4, 4], [-4, 4]],\n",
    "              cmap=\"Blues\")\n",
    "axs[0].set_xlim(-4, 4)\n",
    "axs[0].set_ylim(-4, 4)\n",
    "axs[0].set_xlabel(\"x\")\n",
    "axs[0].set_ylabel(\"p\")\n",
    "axs[0].set_xticks([-4, -2, 0, 2, 4])\n",
    "axs[0].set_yticks([-4, -2, 0, 2, 4])\n",
    "axs[0].set_aspect(1)\n",
    "\n",
    "position = np.linspace(-4, 4, 500)\n",
    "boltz_factor = np.exp(-0.5/lan.kT*lan.k_spring*position**2)\n",
    "Z = np.trapz(boltz_factor, position)\n",
    "boltz_factor /= Z\n",
    "axs[1].plot(position, boltz_factor, color=\"k\", linestyle=\"--\")\n",
    "axs[2].plot(position, boltz_factor, color=\"k\", linestyle=\"--\")\n",
    "\n",
    "axs[1].hist(lan.x_list, density=True, bins=21, alpha=0.8)\n",
    "axs[1].set_xlim(-4, 4)\n",
    "axs[1].set_ylim(0, 0.5)\n",
    "axs[1].set_box_aspect(1)\n",
    "axs[1].set_xlabel(\"x\")\n",
    "axs[1].set_ylabel(\"probability\")\n",
    "axs[1].set_xticks([-4, -2, 0, 2, 4])\n",
    "\n",
    "axs[2].hist(lan.v_list, density=True, bins=21, alpha=0.8)\n",
    "axs[2].set_xlim(-4, 4)\n",
    "axs[2].set_ylim(0, 0.5) \n",
    "axs[2].set_box_aspect(1)\n",
    "axs[2].set_xlabel(\"p\")\n",
    "axs[2].set_ylabel(\"probability\")\n",
    "axs[2].set_xticks([-4, -2, 0, 2, 4])\n",
    "\n",
    "fig.set_size_inches(11, 3)\n",
    "plt.savefig(\"md_lan_nh.pdf\", bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ovito",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
